---
title: "Growth Cone Data"
author: "Adam Bindas"
date: "7/05/2022"
output: html_document
---

```{r setup, include=FALSE}
# Load packages
library(ggplot2)
library(moments)
library(emmeans)
library(lme4)
library(lmtest)
library(ggsignif)

# Import Data (csv from excel)
g_cone <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Growth Cone Live Imaging/Growth cone quantification 22-05-27.csv", header=TRUE, stringsAsFactors=FALSE)
```

```{r Inspection/Selection}
# Define variables for analysis (e.g., group factor, dosage level)
g_cone$Groups.f <- factor(g_cone$group)
g_cone$Groups.f <- factor(g_cone$Groups.f,levels = c("control", "PFF", "4xPFF"))
g_cone$control <- ifelse(g_cone$group == "control",1,0)
g_cone$pff1 <- ifelse(g_cone$group == "PFF",1,0)
g_cone$pff4 <- ifelse(g_cone$group == "4xPFF",1,0)
g_cone$dose <- ifelse(g_cone$group == "PFF",1,0)
g_cone$dose <- ifelse(g_cone$group == "4xPFF",4,g_cone$dose)
g_cone$d21 <- ifelse(g_cone$day == 21,1,0)
g_cone$att <- factor(g_cone$rat)
g_cone$qualz <- factor(g_cone$qual.focal)
g_cone$day.no <- factor(g_cone$day)

# Define Complex Variables that include multiple outcome variables
g_cone$eventpcount <- g_cone$no.events.permin/(g_cone$ext.count)
g_cone$ext.ev.area <- g_cone$no.events.permin/(g_cone$ext.count*g_cone$area.micromet)
g_cone$ext.p.area <- g_cone$no.events.permin/g_cone$area.micromet

# Quality Removal Subset. Quality is defined manually during data acquisition (0 = clear throughout recording, 1 = partially clear but quantifiable, 2 = poor clarity, in cases an estimation may have been conducted)
g_cone2 <- subset(g_cone, qualz != 2)
g_cone14 <- subset(g_cone2, day.no == 14)
g_cone21 <- subset(g_cone2, day.no == 21)
```


``` {r Sample Size}
# Identify sample size that is analyzed. N = experimental replicate, m = wells, d = number of images
# length(unique(g_cone$rat))

# For each level of information, identify the number of unique values (e.g., attempt number) across each day and experimental group
n.count.area <- aggregate(data = g_cone2,
                          att ~ day.no + Groups.f,
                          function(att) length(unique(att)))
m.count.area <- aggregate(data = g_cone2,
                          m.id ~ day.no + Groups.f,
                          function(m.id) length(unique(m.id)))
d.count.area  <- aggregate(data = g_cone2,
                          unique.id ~ day.no + Groups.f,
                          function(unique.id) length(unique(unique.id)))

# Merge output sample size variables into a table form
gcone.replicates <- cbind(n.count.area[,3],m.count.area[,3],d.count.area[,3])

# Label column names
colnames(gcone.replicates) <- c("N","m","no.images")

# Label row names with a combined variable
rownames(gcone.replicates) <- paste(m.count.area[,2]," Day ",m.count.area[,1])

# Display sample size
gcone.replicates
```


```{r Growth Cone Area graphing/analysis}
#Histogram of Growth Cone Area pre- and post-transformation with skewness
hist(g_cone$area.micromet)
hist(g_cone$t.area.micromet <- log(g_cone$area.micromet))
skewness(g_cone$t.area.micromet, na.rm = TRUE)

#Plot Values with each day included
plot.area1 <- ggplot(g_cone, aes(x=Groups.f, y=t.area.micromet, fill=day.no)) +
  ylab("Log of Growth Cone Area (µm2)") + 
  xlab("Groups") + 
  labs(
    title = "Growth Cone Area is not appreciably changed by PFF dosage", 
    subtitle = "N = 4-5, m = 9-12, # cones = 36-54, error bars = SEM", 
    fill = "Replicate")+
  scale_x_discrete(labels=c("Control", "PFF", "4X PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = g_cone, aes(x = Groups.f, y = area.micromet, 
                                  group = interaction(Groups.f, day.no)),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.3,
               stackdir = "center", stackratio = 1)

plot.area1

# Analysis
# To determine whether adjusting for random effects would be valuable, a likelihood ratio test is conducted between models with and without random intercept effects. The simplest test (if any) that is significantly improved from a model without adjustments.

# Day 14 analysis.
m.area14 <- lm(log(area.micromet) ~ 1 + Groups.f, 
                  data = g_cone14)
m.area14.2 <- lmer(log(area.micromet) ~ 1 + Groups.f +
                    (1|att),  REML=FALSE, 
                  data = g_cone14)

#Likelihood ratio test
lrtest(m.area14,m.area14.2)

#Summarize whether there are significant differences between groups
summary(m.area14)

# Summarize model using emmeans for estimated means 
summary(emmeans(m.area14,~Groups.f,type = "response"))
       
# Day 14 analysis.
m.area21 <- lm(log(area.micromet) ~ 1 + Groups.f, 
                  data = g_cone21)
m.area21.2 <- lmer(log(area.micromet) ~ 1 + Groups.f +
                    (1|att),  REML=FALSE, 
                  data = g_cone21)

#Likelihood ratio test
lrtest(m.area21,m.area21.2)

#Summarize whether there are significant differences between groups
summary(m.area21.2)
summary(emmeans(m.area21,~Groups.f,type = "response"))
```

```{r extensions}
#Histogram
hist(g_cone2$ext.count)
hist(g_cone2$t.ext.count <- log(g_cone2$ext.count))
skewness(g_cone2$t.ext.count, na.rm = TRUE)

#Plot Values
plot.ext.count1 <- ggplot(g_cone2, aes(x=Groups.f, y=ext.count, fill=day.no)) +
  ylab("Log of Growth Cone Extension Count") + 
  xlab("Groups") + 
  labs(
    title = "Growth Cone Extention Count is adjusted by PFF dosage in select pairs", 
    subtitle = "N = 4-5, m = 9-12, # cones = 36-54, error bars = SEM", 
    fill = "Replicate")+
  scale_x_discrete(labels=c("Control", "PFF", "4X PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = g_cone2, aes(x = Groups.f, y = ext.count, 
                                  group = interaction(Groups.f, day.no)),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
               stackdir = "center", stackratio = 1)

plot.ext.count1

#Initial Analysis
m.ext.count14 <- lmer(log(ext.count) ~ 1 + Groups.f +
                    (1|att),  REML=FALSE, 
                  data = g_cone14)
m.ext.count14.2 <- lm(log(ext.count) ~ 1 + Groups.f, 
                  data = g_cone14)
anova(m.ext.count14,m.ext.count14.2)
summary(m.ext.count14.2)
summary(em.ext.count14 <- emmeans(m.ext.count14.2, ~Groups.f,type = "response"))

m.ext.count21.1 <- lmer(log(ext.count) ~ 1 + Groups.f +
                    (1|att), REML=FALSE, 
                  data = g_cone21)
m.ext.count21.2 <- lm(log(ext.count) ~ 1 + Groups.f,
                  data = g_cone21)
anova(m.ext.count21.1,m.ext.count21.2)
summary(m.ext.count21.2)
summary(em.ext.count21 <- emmeans(m.ext.count21.2, ~Groups.f,type = "response"))
```


```{r event rate}
#Histogram
hist(g_cone2$no.events.permin)
hist(g_cone2$t.events <- (g_cone2$no.events.permin)^(1/4))
skewness(g_cone2$t.events, na.rm = TRUE)

#Plot Values
plot.events1 <- ggplot(g_cone2, aes(x=Groups.f, y=no.events.permin, fill=day.no)) +
  ylab("Fourth Root of Number of Events per Image") + 
  xlab("Groups") + 
  labs(
    title = "Growth Cone Event rate is trending towards minimally being impacted by PFF dosage", 
    subtitle = "N = 4-5, m = 9-12, # cones = 36-54, error bars = SEM", 
    fill = "Replicate")+
  scale_x_discrete(labels=c("Control", "PFF", "4X PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = g_cone2, aes(x = Groups.f, y = no.events.permin, 
                                  group = interaction(Groups.f, day.no)),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
               stackdir = "center", stackratio = 1)

plot.events1

m.event.rate14 <- lmer(no.events.permin^0.25 ~ 1 + Groups.f +
                    (1|att), 
                  data = g_cone14)
m.event.rate14.2 <- lm(no.events.permin^0.25 ~ 1 + Groups.f, 
                  data = g_cone14)
m.event.rate14.3 <- lmer(no.events.permin^0.25 ~ 1 + Groups.f +
                    (1|att) + (1|att:well),
                  data = g_cone14)
anova(m.event.rate14,m.event.rate14.2,m.event.rate14.3)
summary(m.event.rate14.2)

event.rate14.tran <- make.tran("power", 0.25)
warp.event.rate14 <- with(event.rate14.tran, 
    lm(linkfun(no.events.permin) ~ Groups.f, data = g_cone14))
summary(em.event.rate14 <- emmeans(warp.event.rate14, ~Groups.f, type = "response"))

m.event.rate21 <- lm(no.events.permin^0.25 ~ 1 + Groups.f, 
                  data = g_cone21)
m.event.rate21.2 <- lmer(no.events.permin^0.25 ~ 1 + Groups.f +
                    (1|att), 
                  data = g_cone21)
m.event.rate21.3 <- lmer(no.events.permin^0.25 ~ 1 + Groups.f +
                    (1|att) + (1|att:well),
                  data = g_cone21)
lrtest(m.event.rate21,m.event.rate21.2,m.event.rate21.3)
summary(m.event.rate21.2)
warp.event.rate21 <- with(event.rate14.tran, 
    lmer(linkfun(no.events.permin) ~ Groups.f + (1|att), data = g_cone21))
summary(em.event.rate21 <- emmeans(warp.event.rate21, ~Groups.f, type = "response"))
```

```{r movements per slide and extention}
#Histogram
hist(g_cone2$eventpcount)
hist(g_cone2$t.eventpcount <- (g_cone2$eventpcount)^(1/4))
skewness(g_cone2$t.eventpcount, na.rm = TRUE)

#Plot Values
plot.eventpcount1 <- ggplot(g_cone2, aes(x=Groups.f, y=eventpcount, fill=day.no)) +
  ylab("Square of root event rate per extension") + 
  xlab("Groups") + 
  labs(
    title = "Growth Cone Event rate per extension is trending towards minimally being impacted by PFF dosage", 
    subtitle = "N = 4-5, m = 9-12, # cones = 36-54, error bars = SEM", 
    fill = "Replicate")+
  scale_x_discrete(labels=c("Control", "PFF", "4X PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = g_cone2, aes(x = Groups.f, y = eventpcount, 
                                  group = interaction(Groups.f, day.no)),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
               stackdir = "center", stackratio = 1)

plot.eventpcount1

m.rate.p.count14 <- lm(eventpcount^0.25 ~ 1 + Groups.f,
                  data = g_cone14)
m.rate.p.count14.2 <- lmer(eventpcount^0.25 ~ 1 + Groups.f + 
                    (1|att), 
                  data = g_cone14)
m.rate.p.count14.3 <- lmer(eventpcount^0.25 ~ 1 + Groups.f + 
                    (1|att) + (1|att:well), 
                  data = g_cone14)
lrtest(m.rate.p.count14,m.rate.p.count14.2,m.rate.p.count14.3)
lrtest(m.rate.p.count14.2,m.rate.p.count14.3)
summary(m.rate.p.count14.3)

m.rate.p.count21 <- lm(eventpcount^0.25 ~ 1 + Groups.f,
                  data = g_cone21)
m.rate.p.count21.2 <- lmer(eventpcount^0.25 ~ 1 + Groups.f + 
                    (1|att), 
                  data = g_cone21)
m.rate.p.count21.3 <- lmer(eventpcount^0.25 ~ 1 + Groups.f + 
                    (1|att) + (1|att:well), 
                  data = g_cone21)
lrtest(m.rate.p.count21,m.rate.p.count21.2,m.rate.p.count21.3)
summary(m.rate.p.count21.2)
event.rate21.tran <- make.tran("power", 0.25)
warp.event.rate.p.count21 <- with(event.rate21.tran, 
    lmer(linkfun(eventpcount) ~ Groups.f + (1|att), data = g_cone21))
summary(em.event.rate.p.count21 <- emmeans(warp.event.rate.p.count21, ~Groups.f, type = "response"))
```
